#include "SmithWaterman.h"

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
#define max(x, y)        (x)>(y)?(x):(y)
#define max1(x, y, z)    max(max(x,y), z)
#define max2(x, y, z, k) max(max(max((x),(y)), (z)), (k))

#define USEMATCHTAG            1
#define GAP_CHARACTER         '-'
#define GAP_TAG               ' '
#define MISMATCH_TAG          ' '
#define APPROXIMATE_MATCH_TAG '+'
#define MATCH_TAG             '|'
#define OUTPUT                 1
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
namespace MatrixAlignmentSpace {

  bool FindMaxScoreCell(TNode <MaxScoreCell>* NodeT, void* ref = NULL) {
    MaxScoreCell *compareCell = (TNode <MaxScoreCell>*)ref;
    MaxScoreCell *listCell    = NodeT->Info();
    return ((compareCell->maxIScore == listCell->maxIScore) && (compareCell->maxIRow == listCell->maxIRow) && (compareCell->maxICol == listCell->maxICol))?true:false;
  }

};
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
SmithWaterman::SmithWaterman() {
    
  this->seq1             = new String();
  this->seq2             = new String();
  this->alignment        = new String();

  this->regionsConserved = new TListE <String> ();

  this->numberRegions    = 0;
  this->matrix           = NULL;

  this->max_row          = 0;
  this->max_col          = 0;
  this->score            = 0;

  this->listMaxCell      = new TListE <String> ();
  
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
SmithWaterman::SmithWaterman(String* seq1, String* seq2):score(0), numberRegions(0) {

  this->seq1 = new String();
  this->seq2 = new String();

  if (seq1) this->seq1->In(seq1);
  if (seq2) this->seq2->In(seq2);
  
  this->max_row = this->seq1->Length();
  this->max_col = this->seq2->Length();
  
  this->matrix = new int*[this->max_row + 1];
  for (int i = 0; i <= this->max_row; i++) {
    this->matrix[i] = new int[this->max_col + 1];
    for (int j = 0; j <= this->max_col; j++)
      this->matrix[i][j] = 0;
  }

  this->alignment        = new String();
  this->regionsConserved = new TListE <String> ();
  
  this->numberRegions    = 0;
  this->score            = 0;

  this->listMaxCell      = new TListE <String> ();
  
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
SmithWaterman::SmithWaterman(const SmithWaterman& smithWaterman) {

  this->seq1 = new String();
  this->seq2 = new String();

  if (smithWaterman.seq1->Length()) this->seq1->In(smithWaterman.seq1);
  if (smithWaterman.seq2->Length()) this->seq2->In(smithWaterman.seq2);

  this->max_row = this->seq1->Length();
  this->max_col = this->seq2->Length();

  this->matrix = new int*[this->max_row + 1];
  for (int i = 0; i <= this->max_row; i++) {
    this->matrix[i] = new int[this->max_col + 1];
    for (int j = 0; j <= this->max_col; j++)
      this->matrix[i][j] = 0;
  }
  
  this->score     = smithWaterman.score;

  this->alignment = new String(smithWaterman.aligment);

  if (alignmentT.alignment) {
    alignment = new char[strlen(alignmentT.alignment) + 1];
    strcpy(alignment, alignmentT.alignment);
  }
  else alignment = NULL;

  if (alignmentT.numberRegions && alignmentT.regionsConserved) {
      
    regionsConserved = new char*[alignmentT.numberRegions];
    
    for (int i = 0; i < alignmentT.numberRegions; i++) {
      if (alignmentT.regionsConserved[i]) {
        regionsConserved[i] = new char[strlen(alignmentT.regionsConserved[i]) + 1];
        strcpy(regionsConserved[i], alignmentT.regionsConserved[i]);
      }
      else regionsConserved[i] = NULL;
    }
    
  }
  else regionsConserved = NULL;
  
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
SmithWaterman::~SmithWaterman() {
    
  if (seq1) delete[] seq1;
  if (seq2) delete[] seq2;

  for (int i = 0; i <= max_row; i++)
    if (matrix[i]) delete[] matrix[i];
  if (matrix) delete[] matrix;

  if (alignment) delete[] alignment;

  for (int i = 0; i < numberRegions; i++)
    if (regionsConserved[i]) delete[] regionsConserved[i];
  if (regionsConserved) delete[] regionsConserved;

  if (listMaxCell) delete listMaxCell;
  
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
void SmithWaterman::SequencesIn(const char* seq1In, const char* seq2In) {

  this->CleanSequences();
    
  if (seq1In) {
    seq1 = new char[strlen(seq1In) + 1];
    strcpy(seq1, seq1In);
    seq1[strlen(seq1In)] = '\0';
  }
  else seq1 = NULL;
  if (seq2In) {
    seq2 = new char[strlen(seq2In) + 1];
    strcpy(seq2, seq2In);
    seq2[strlen(seq2In)] = '\0';
  }
  else seq2 = NULL;

  max_row = (seq1)?strlen(seq1):0;
  max_col = (seq2)?strlen(seq2):0;

  matrix = new int*[max_row + 1];
  for (int i = 0; i <= max_row; i++) {
    matrix[i] = new int[max_col + 1];
    for (int j = 0; j <= max_col; j++)
      matrix[i][j] = 0;
  }
  
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
void SmithWaterman::SequencesOut(char* seq1In, char* seq2In) {
    
  if (seq1In) delete[] seq1In;
  if (seq2In) delete[] seq2In;

  if (seq1) {
    seq1In = new char[strlen(seq1) + 1];
    strcpy(seq1In, seq1);
  }
  else seq1In = NULL;
  if (seq2) {
    seq2In = new char[strlen(seq2) + 1];
    strcpy(seq2In, seq2);
  }
  else seq2In = NULL;

}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
float SmithWaterman::Score() {
  return score;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
void SmithWaterman::CleanSequences(void){
    
  if (seq1)      delete[] seq1;
  if (seq2)      delete[] seq2;
  if (alignment) delete[] seq2;

  if (max_row)
    for (int i = 0; i <= max_row; i++)
      if (matrix[i]) delete[] matrix[i];
    
  if (matrix) delete[] matrix;

  for (int i = 0; i < numberRegions; i++)
    if (regionsConserved[i]) delete[] regionsConserved[i];
  delete[] regionsConserved;

  max_row       = 0;
  max_col       = 0;
  score         = 0;
  numberRegions = 0;
  
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
char* SmithWaterman::ComputePairwiseAlignment (){

  this->ComputeMatrix();
  alignment = this->BuildOptimalAlignment();
  return alignment;
  
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
void SmithWaterman::ComputeMatrix() {
    
  int rows, cols, ins, sub, del, max_score;

  rows = strlen(seq1) + 1;
  cols = strlen(seq2) + 1;

  // keep track of the maximum score
   this->max_row = this->max_col = max_score = 0;

  // calculates the similarity matrix (row-wise)
  for (int i = 1; i < rows; i++) {
      
    // initiate first column
    matrix[i][0] = 0;

    for (int j = 1; j < cols; j++) {
        
      ins = matrix[i    ][j - 1] + this->ScoreInsertion(seq2[j - 1]);
      sub = matrix[i - 1][j - 1] + this->ScoreSubstitution(seq1[i - 1], seq2[j - 1]);
      del = matrix[i - 1][j    ] + this->ScoreDeletion(seq1[i - 1]);

       // choose the greatest
       matrix[i][j] = max2(ins, sub, del, 0);
       if (matrix[i][j] > max_score) {
          // keep track of the maximum score
          max_score     = matrix[i][j];
          this->max_row = i;
          this->max_col = j;
        }
       
      }
    }

  if (OUTPUT) {
    printf("-\t\t-");
    for (int i = 0; i < max_col; i++) printf("\t\t%c", seq2[i]);
    printf("\n");
    for (int i = 0; i < rows; i++) {
      if (i) printf("%c\t\t", seq1[i - 1]);
      else   printf("-\t\t");
      for (int j = 0; j < cols; j++)
        printf("[%d, %d]=%d\t", i, j, matrix[i][j]);
      printf("\t\n");
    }
  }

}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
char* SmithWaterman::BuildOptimalAlignment(void) {

  char *gapped_seq1 = NULL, *score_tag_line = NULL, *gapped_seq2 = NULL;
  int	i, j, max_score, sub;

  //start at the cell with maximum score
  i = this->max_row;
  j = this->max_col;

  max_score = matrix[i][j];

  printf("Optimal Aligment Begin in [%d, %d]\n", i, j);

  gapped_seq1    = new char[strlen(seq1) + strlen(seq2) + 1];
  score_tag_line = new char[strlen(seq1) + strlen(seq2) + 1];
  gapped_seq2	 = new char[strlen(seq1) + strlen(seq2) + 1];

  int gapped_seq1C = 0, score_tag_lineC = 0, gapped_seq2C = 0;

  gapped_seq1[gapped_seq1C]         = '\0';
  score_tag_line[score_tag_lineC]   = '\0';
  gapped_seq2[gapped_seq2C]         = '\0';
  
  while ((i > 0) || (j > 0) && (matrix[i][j] > 0)) {
        
    if (j > 0)
      if (matrix[i][j] == matrix[i][j - 1] + this->ScoreInsertion(seq2[j - 1])) {
            
        //insertion
        gapped_seq1[gapped_seq1C++]       = GAP_CHARACTER;
        score_tag_line[score_tag_lineC++] = GAP_TAG;
        gapped_seq2[gapped_seq2C++]       = seq2[j - 1];

        j = j - 1;

        //skip to the next iteration
        continue;
      }

    if ((i > 0) && (j > 0)) {
          
      sub = this->ScoreSubstitution(seq1[i - 1], seq2[j - 1]);

      if (matrix[i][j] == matrix[i - 1][j - 1] + sub) {
            
        //substitution
        gapped_seq1[gapped_seq1C++] = seq1[i - 1];
        if (seq1[i - 1] == seq2[j - 1])
          if   (USEMATCHTAG) score_tag_line[score_tag_lineC++] = MATCH_TAG;
          else               score_tag_line[score_tag_lineC++] = seq1[i - 1];
        else if (sub > 0)
          score_tag_line[score_tag_lineC++] = APPROXIMATE_MATCH_TAG;
        else
          score_tag_line[score_tag_lineC++] = MISMATCH_TAG;
          
        gapped_seq2[gapped_seq2C++] = seq2[j - 1];

        i--; j--;

        // skip to the next iteration
        continue;
      }
    }

    // must be a deletion
    gapped_seq1[gapped_seq1C++]       = seq1[i - 1];
    score_tag_line[score_tag_lineC++] = GAP_TAG;
    gapped_seq2[gapped_seq2C++]       = GAP_CHARACTER;

    i = i - 1;
    
  }
  
  gapped_seq1[gapped_seq1C]       = '\0';
  score_tag_line[score_tag_lineC] = '\0';
  gapped_seq2[gapped_seq2C]       = '\0';

  InvertSeq(gapped_seq1); InvertSeq(score_tag_line); InvertSeq(gapped_seq2);
  
  char *result = new char[strlen(gapped_seq1) + strlen(score_tag_line) + strlen(gapped_seq2) + 4];
  sprintf(result, "%s\n%s\n%s\0", gapped_seq1, score_tag_line, gapped_seq2);

  delete[] gapped_seq1;
  delete[] score_tag_line;
  delete[] gapped_seq2;

  return result;
  
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int SmithWaterman::ComputeScore(void) {
    
  int	*array;
  int 	 rows = strlen(seq1) + 1, cols = strlen(seq2) + 1;
  int 	 i, j, tmp, ins, del, sub, max_score;

  //keep track of the maximum score
  max_score = 0;

  if (rows <= cols) {
      
    // goes columnwise
    array = new int[rows];

    // initiate first column
    for (i = 0; i < rows; i++) array[i] = 0;

    // calculate the similarity matrix (keep current column only)
    for (j = 1; j < cols; j++) {

      // set first position to zero (tmp hold values
      // that will be later moved to the array)
      tmp = 0;

      for (i = 1; i < rows; i++) {
          
        ins = array[i    ] + this->ScoreInsertion(seq2[j - 1]);
        sub = array[i - 1] + this->ScoreSubstitution(seq1[i - 1], seq2[j - 1]);
        del = tmp          + this->ScoreDeletion(seq1[i - 1]);

        // move the temp value to the array
        array[i - 1] = tmp;

        // choose the greatest (or zero if all negative)
        tmp = max2(ins, sub, del, 0);

        // keep track of the maximum score
        if (tmp > max_score) max_score = tmp;
      }

      // move the temp value to the array
      array[rows - 1] = tmp;
    }
  }
  else {
      
    // goes rowwise
    array = new int [cols];

    // initiate first row
    for (j = 0; j < cols; j++) array[j] = 0;

    // calculate the similarity matrix (keep current row only)
    for (i = 1; i < rows; i++) {
        
      // set first position to zero (tmp hold values
      // that will be later moved to the array)
      tmp = 0;

      for (j = 1; j < cols; j++) {
          
        ins = tmp          + this->ScoreInsertion(seq2[j - 1]);
        sub = array[j - 1] + this->ScoreSubstitution(seq1[i - 1], seq2[j - 1]);
        del = array[j    ] + this->ScoreDeletion(seq1[i]);

        // move the temp value to the array
        array[j - 1] = tmp;

        // choose the greatest (or zero if all negative)
        tmp = max2(ins, sub, del, 0);

        // keep track of the maximum score
        if (tmp > max_score) max_score = tmp;
        
      }

      // move the temp value to the array
      array[cols - 1] = tmp;
    }
  }

  delete[] array;

  return max_score;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int SmithWaterman::ScoreInsertion(char elem) {
  return -3;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int SmithWaterman::ScoreSubstitution(char elem1, char elem2) {
  return (elem1 == elem2)?10:-10;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
int SmithWaterman::ScoreDeletion(char elem) {
  return -3;
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
void SmithWaterman::InvertSeq(char *Sequence) {

  for (int i = 0; i < strlen(Sequence)/2; i++) {
    char character                     = Sequence[i];
    Sequence[i]                        = Sequence[strlen(Sequence) - i - 1];
    Sequence[strlen(Sequence) - i - 1] = character;
  }    

}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
char* SmithWaterman::FindBestStrictLocalAlignment(void) {

  this->ComputeMatrix();
  this->FindMaxScoreIterative();

  char *gapped_seq1 = NULL, *score_tag_line = NULL, *gapped_seq2 = NULL;
  int	i, j, max_score, sub;

  //start at the cell with maximum score
  i = this->max_row;
  j = this->max_col;

  max_score = matrix[i][j];

  gapped_seq1    = new char[strlen(seq1) + strlen(seq2) + 1];
  score_tag_line = new char[strlen(seq1) + strlen(seq2) + 1];
  gapped_seq2	 = new char[strlen(seq1) + strlen(seq2) + 1];

  int gapped_seq1C = 0, score_tag_lineC = 0, gapped_seq2C = 0;

  gapped_seq1[gapped_seq1C]         = '\0';
  score_tag_line[score_tag_lineC]   = '\0';
  gapped_seq2[gapped_seq2C]         = '\0';

  while ((j > 0 || j > 0) && (matrix[i][j] > 0)) {

    if (j > 0)
      if (matrix[i][j] == matrix[i][j - 1] + this->ScoreInsertion(seq2[j - 1])) {

        //insertion
        gapped_seq1[gapped_seq1C++]       = GAP_CHARACTER;
        score_tag_line[score_tag_lineC++] = GAP_TAG;
        gapped_seq2[gapped_seq2C++]       = seq2[j - 1];

        j = j - 1;

        //skip to the next iteration
        continue;
      }

    if ((i > 0) && (j > 0)) {

      sub = this->ScoreSubstitution(seq1[i - 1], seq2[j - 1]);

      if (matrix[i][j] == matrix[i - 1][j - 1] + sub) {

        //substitution
        gapped_seq1[gapped_seq1C++] = seq1[j - 1];
        if (seq1[i - 1] == seq2[j - 1])
          if   (USEMATCHTAG) score_tag_line[score_tag_lineC++] = MATCH_TAG;
          else               score_tag_line[score_tag_lineC++] = seq1[i - 1];
        else if (sub > 0)
          score_tag_line[score_tag_lineC++] = APPROXIMATE_MATCH_TAG;
        else
          score_tag_line[score_tag_lineC++] = MISMATCH_TAG;

        gapped_seq2[gapped_seq2C++] = seq2[j - 1];

        i = i - 1; j = j - 1;

        // skip to the next iteration
        continue;
      }
    }

    // must be a deletion
    gapped_seq1[gapped_seq1C++]       = seq1[i - 1];
    score_tag_line[score_tag_lineC++] = GAP_TAG;
    gapped_seq2[gapped_seq2C++]       = GAP_CHARACTER;

    i = i - 1;

  }

  gapped_seq1[gapped_seq1C]       = '\0';
  score_tag_line[score_tag_lineC] = '\0';
  gapped_seq2[gapped_seq2C]       = '\0';

  InvertSeq(gapped_seq1); InvertSeq(score_tag_line); InvertSeq(gapped_seq2);

  char *result = new char[strlen(gapped_seq1) + strlen(score_tag_line) + strlen(gapped_seq2) + 4];
  sprintf(result, "%s\n%s\n%s\0", gapped_seq1, score_tag_line, gapped_seq2);

  delete[] gapped_seq1;
  delete[] score_tag_line;
  delete[] gapped_seq2;

  return result;

}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
void SmithWaterman::FindMaxScoreIterative(void) {

  int i = 0, scoreBase = 0, posX = 0, posY = 0;

  if (listMaxCell) delete listMaxCell;
  listMaxCell = new TListE;

  do {

    int maximunScore = 0, row = 0, col = 0;
    int initialRow = max_row;
    int initialCol = max_col;

    for (int i = initialRow; i > 0; i--)
      for (int j = initialCol; j > 0; j--) {
        if ((matrix[i][j] > maximunScore) && (seq1[i - 1] == seq2[j - 1])) {
          if (scoreBase > 0) {

            MaxScoreCell compareCell;
            
            compareCell.maxIScore = matrix[i][j];
            compareCell.maxIRow   = i;
            compareCell.maxICol   = j;

            MaxScoreCell *tmpMaxCell = (MaxScoreCell*)(listMaxCell->FirstOcurrence(MatrixAlignmentSpace::FindMaxScoreCell, &compareCell));

            if ((i == posX) && (j == posY)) continue;
            else if (tmpMaxCell) continue;
            else if (matrix[i][j] <= scoreBase){
              row           = i;
              col           = j;
              maximunScore  = matrix[i][j];
            }
          }
          else {
            row           = i;
            col           = j;
            maximunScore  = matrix[i][j];
          }
        }
      }

    MaxScoreCell *tmpMaxCell = new MaxScoreCell;
    
    posX      = tmpMaxCell->maxIRow   = row;
    posY      = tmpMaxCell->maxICol   = col;
    scoreBase = tmpMaxCell->maxIScore = maximunScore;

    if (scoreBase) listMaxCell->Add(tmpMaxCell);
    
  } while (scoreBase);

    
}
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
//---------------------------------------------------------------------------------------------------------------------------------------------------------//
